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Abstract 

We investigate the time evolution of a classical ensemble of isolated peri- 
odic chains of 0(iV)-symmetric anharmonic oscillators. Our method is based 
on an exact evolution equation for the time dependence of correlation func- 
tions. We discuss its solutions in an approximation which retains all con- 
tributions in next-to- leading order in a 1/N expansion and preserves time 
reflection symmetry. We observe effective irreversibility and approximate 
thermalization. At large time the system approaches stationary solutions 
in the vicinity of, but not identical to, thermal equilibrium. The ensemble 
therefore retains some memory of the initial condition beyond the conserved 
total energy. Such a behavior with incomplete thermalization is referred to 
as "mesoscopic dynamics". It is expected for systems in a small volume. 
Surprisingly, we find that the nonthermal asymptotic stationary solutions do 
not change for large volume. This raises questions on Boltzmann's conjecture 
that macroscopic isolated systems thermalize. 
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1 Introduction 



A central piece in our understanding of the dynamics of large statistical 
systems is Boltzmann's conjecture that an ensemble of isolated interacting 
systems approaches thermal equilibrium at large times. The asymptotic val- 
ues of the correlation functions can then be computed from the resulting 
microcanonical equilibrium ensemble. According to Boltzmann's conjecture 
their values only depend on the energy density of the system or, equivalently, 
the temperature (with standard modifications in the presence of other con- 
served extensive quantities). Besides the energy, all memory about the initial 
conditions is lost asymptotically - and also in practice if typical relaxation 
time scales are not too large. The thermalization conjecture only applies to 
spatially extended systems in the limit of infinite volume. So far no proof 
for this hypothesis has been given. Correspondingly, the question of how 
effective irreversibility arises from microscopic equations which are invariant 
under time reflection, or from the time reversible Liouville equation, has not 
found a complete answer todate. This effective irreversibility is the basis of 
widely used effective equations (like the Boltzmann equation). 

We want to address two issues by a direct study of the time dependence 
of the correlation functions. The first concerns isolated systems with a finite 
number of degrees of freedom, corresponding to a finite volume V. No equili- 
bration is expected for microscopic systems of only a few degrees of freedom. 
Recently, this has been demonstrated explicitly for the correlation functions 
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of coupled anharmonic oscillators PJ []. Let us assume for a moment that 
thermalization occurs for macroscopic systems. Then a smooth transition 
between the two extremes requires that there must be some intermediate 
volume or number of degrees of freedom where thermalization is incomplete, 
i.e. the ensemble retains some memory of the initial conditions. We may 
call the associated time evolution of the correlation functions "mesoscopic 
dynamics". Mesoscopic dynamics is characterized by partial effective irre- 
versibility on one side, i.e. many details of the initial conditions get lost for 
asymptotic times. On the other hand, this "loss of memory" is not complete 
(as for strictly thermalizing systems) so that partial information about the 
initial state can be recovered even after an arbitrarily long time. 

The second question concerns the validity of Boltzmann's conjecture for 
systems with a large volume. It is conceivable that thermalization remains 
incomplete even for macroscopic isolated systems. In this case some charac- 
teristics of mesoscopic dynamics would survive in the infinite volume limit. 
We emphasize that mesoscopic dynamics is always relevant in an appropriate 
volume range. Our second question therefore asks if and how the asymp- 
totic loss of memory becomes complete as the volume becomes macroscopic. 
Within our approximations we find that certain features of mesoscopic dy- 
namics remain present for isolated systems in the large volume limit! 

Our investigation is based on an exact evolution equation for the time 

dependent effective action which is the generating functional of the equal 

lr The issue would be completely different, of course, if the system were coupled to an 
external heat bath. 
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time one particle irreducible (1PI) correlation functions. The direct study of 
the time evolution of the correlation functions circumvents the calculation of 
the time dependence of the probability distribution or density matrix. For 
classical systems the exact evolution equation is equivalent to the BBGKY 
hierarchy Q. A study of the 1PI correlation functions offers, however, new 
possibilities of systematic truncations. In particular, thermal equilibrium is 
now present as a stationary solution at every step of the truncation. Also 
the generalization to quantum statistics is straightforward and surprisingly 
simple 0. An investigation of the general structure of the exact evolution 
equation for the time dependent effective action reveals many new stationary 
solutions besides thermal equilibrium |J. The question of their dynamical 
role is part of the scope of this paper. 

We will concentrate here on a particular example, namely a periodic chain 
of coupled anharmonic oscillators with O(N) symmetry. The Hamiltonian 



H 



dx 



~Pa{x)Pa{x) + ^q a (x)(m 2 - A)q a (x) + (q a (x)q a (x)) 2 



(1) 



describes an interacting time reversible system with microscopic time evolu- 
tion given by 



d t q a {x) = p a (x) 



A 



dtPa(x) = -(m 2 - A)q a (x) - ^q b (x)q b (x)q a (x). (2) 
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We may assume that the oscillators sit on discrete lattice sites with dis- 
tance a such that high momenta are cut off. Then A is an appropriate dis- 
cretized Laplace operator involving neighboring sites. The index a = 1...N 
counts the oscillators at a given site and H preserves O(N) symmetry. Prac- 
tical examples for small N may be found in the form of ring-sized molecules 
with discrete translational symmetry along the ring. For example, the q a 
may describe displacements from the equilibrium position (with p a the asso- 
ciated momenta). The limit of large / also describes large linear molecules if 
boundary effects from the ends can be neglected. Another interesting limit 
is A — > oo for fixed m 2 /X = —1. This imposes the nonlinear constraint 
QaQa — 1 and describes classical bosonic spin chains (N = 3 for spin one). 
The length / of the chain plays the role of the volume in three dimensional 
systems. For I = we are left with an ensemble of simple iV-component 
anharmonic oscillators. This case has been studied in detail in PJ. No ther- 
malization is possible since infinitely many conserved correlation functions 
keep the memory of the initial condition and obstruct thermalization. The 
conserved cumulants are simply related to powers of the conserved energy 
and squared 0(iV)-angular momentum L 2 , i.e. < E r (L 2 ) s > does not change 
in time for arbitrary r and s. If Boltzmann's conjecture is true, asymptotic 
thermalization governs the behavior for I — > oo. In this case we conclude 
that there must be a range of I with mesoscopic dynamics, describing the 
transition between the limits / — > and / — > oo. On the other hand, if 
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Boltzmann's conjecture does not hold for isolated systems of this type some 
features of mesoscopic dynamics are expected to be relevant also for / — ► oo. 

Our system is also a prototype for classical and quantum field theories. 
From this point of view it describes a one-dimensional 0(iV)-symmetric scalar 
field theory. There is no conceptual problem in its generalization to three 
dimensions where it would be relevant for cosmology (e.g. inflation and para- 
metric resonance or dynamical scalar fields playing a role in late cosmology), 
for particle physics (e.g. pions in heavy ion collisions) or for many systems 
in statistical mechanics. 

Although most practical applications are for small N, we also discuss 
the limit of large N. The reason is that one of our truncation schemes is 
a systematic expansion in powers of 1/N. The leading order in the 1/N 
expansion has been discussed by various groups with different methods |J, 
0. In this first approximation infinitely many conserved quantities preclude 
thermalization ||. We include here all contributions in next to leading order 
in 1/N. In particular, this includes scattering in three dimensional field 
theories. Particle numbers for individual momentum modes are no longer 
conserved and there is no immediately visible obstruction to thermalization 
anymore. 

We concentrate in this paper on classical statistics. This has the ad- 
vantage that our results can easily be compared with other methods. In 
particular, it should be feasible to solve the microscopic equations numer- 
ically with given initial values and then take averages over an ensemble of 
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initial values. In this way the equal time correlation functions discussed in 
this paper can be directly measured at any later time. The generalization to 
quantum statistics is straightforward in our approach and will be postponed 
to a subsequent paper. 

We investigate homogeneous (translation invariant) and O(N) symmetric 
ensembles. Individual members of these ensembles do not, of course, share 
this high degree of symmetry. For generic initial conditions the solution of 
the microscopic evolution equation (^) is highly inhomogeneous and shows 
no O(N) symmetry. The high symmetry of the ensemble only means that 
we weight the initial conditions according to a probability distribution that 
exhibits this symmetry. In practice, there is actually no need to specify the 
probability distribution at the initial time to explicitely. It is often more 
effective to specify the correlation functions at to- These will constitute 
the initial data for our differential flow equations. In the present paper we 
mainly consider gaussian initial perturbations from equilibrium, where all 1PI 
n-point functions except the two-point functions are equal to their thermal 
values. 

We find effective irreversibility as a property of the solutions of our time 
reversible flow equation. In a wider sense this is due to the existence of at- 
tractive fixed point solutions. In our context an attractive fixed point does 
not necessarily mean that all solutions for a given class of initial ensembles 
asymptotically reach this fixed point. The characteristic behavior for large t 
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is rather characterized by high-frequency oscillations of the correlation func- 
tions around time averaged stationary mean values. The approach to the 
stationary behavior for the time averaged correlation functions shows three 
characteristic features: 

• First we find a fast initial irreversible behavior on typical microscopic 
time scales between the inverse momentum cutoff A -1 = a/n (with a 
the lattice distance) and the inverse mass mT 1 . An example is given 
in Fig.[l| where we plot the time evolution of the ratio between kinetic 
energy and total energy for an out-of-equilibrium system. Comparison 
with the leading order in the 1/N expansion reveals no qualitative 
difference. We conclude that this first period of irreversibility is not 
related to scattering but rather described by dephasing P, pj. A 
"rough thermalization" takes already place at this very early stage. 

• The first stage of rapid "rough thermalization" does not bring the two- 
point functions near the equilibrium values. In Fig.^| we display the 
evolution of the time-averaged two-point function B(q) which charac- 
terizes the gaussian part of the probability distribution for the momenta 
Pa(q) = Jdx e~ iqx p a (x) by < p a (q)pb(q') > = 2?r% + q')8 ah G™{q), 
G™(q) = - C 2 (q)/(A(q)B(q)))- 1 (see below). In thermal 
equilibrium one expects the Maxwell velocity distribution with B(q) = 
(3 = 1/T independent of q and C(q) = 0. One observes that an initially 
disturbed B(q) approaches a stationary value only on time scales much 
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larger than m _1 . "Scattering" is essential for this aspect of irreversibil- 
ity. This can be seen by a comparison with the leading 1/N behavior 
where B(q) oscillates around a time independent value for every q. In 
leading order 1/N no energy is exchanged between the different Fourier 
modes. This explains why time averaged values for B(q) are stationary 
from the beginning and therefore cannot equilibrate. 

• The exchange of energy between different Fourier modes in next-to- 
leading order in the 1/N expansion drives the time averaged velocity 
distribution toward a stationary value. It may be a surprise that in 
general this stationary value differs from thermal equilibrium with uni- 
form B(q) = p. The implications of this finding will be discussed in 
the last section. 

2 The method 

Our investigation is based on the time-dependent effective action [0, which 
generates the equal-time 1PI correlation functions. We consider an TV-component 
(1 + l)-dimensional scalar 4 theory and ensembles which are invariant un- 
der internal O(N) transformations, spatial translations, and reflection. Our 
truncation retains all 1PI n-point functions up to n = 4 and omits 1PI ver- 
tices with n > 6. One should notice that this still includes connected n-point 
functions with arbitrary n. In this approximation the effective action of our 
model becomes: 
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Figure 1: Energy "equipartition" . We show e^/e as a function of time, in the 
leading (upper panel) and next-to-leading 1/N approximations, for A = 2, 
iV = 20, A = 5, and I = 20.1. The system is initially displaced from 
equilibrium according to a gaussian perturbation in B(q) (/3 = 0.5, Db = 
0.25, Ab = 0.469, qs = 2.5; cf. Eq.flTBD). Horizontal lines correspond to 
thermal equilibrium for T = 2. 
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Figure 2: Time evolution of the two-point correlation function. We plot 
B(A/2) for the same system as in Figfj], in the leading (upper curve) and 
next-to- leading 1/N approximations. The plotted values are averages over 
the time interval [t — 20, t]. Note that the equilibrium value for the corre- 
sponding energy is B(A/2) = 0.5. In the leading 1/N approximation B(q) 
stays esssentially constant for arbitrarily long time. 
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r[0,M = \ I ^{A{qW a {q)Uq) + B{q)K{q>a{q) 



2 J (2n) D 

\u(qi, 02, q-z)<Pa{qi)<Pa{q2)<i>b{q?)<Pb{-qi -02- 
+v(qi, 02, g3)vra(gi)0a(g2)06(g3)0fe(-gi - 02 - gs) 
+w(gi, g 2 , g3)7Ta(gi)7Ta(g2)06(g3)0; ) (-gi - ?a - 93) 
+s(?i, Q2, qs) [K a {.qi)'K b (q2)<Pa{.q-z)<Pb{.-qi -02- qi) 

-^a(0lW(02)<Pb(ls)(pb(-<ll ~ 92 - 93)] (3) 

+y(?i, 92, g3)7ra(gi)vra(g2)vr6(g 3 )0fe(-9i - 02 - 03) 
+z(qi, 02, Oi)^a{qi)^a{q2)^b{q^)'r(b{-qi ~q%- qa)}, 

where the "couplings" A(g), w(gi, 02,0z)i etc. depend on time. The 1PI n- 
point functions are obtained by taking derivatives of V with respect to 
and 7r, the second derivative being the inverse propagator. For example, the 
connected two-point function for q a reads 

< 0a(x)0b(y) >c= G(x - y)5 ab = [ S^n G ^y q{X ~ V) ^ (4) 



where 



(27T) 



G{q) ~ A(q)B(q)-(P( q y (5) 



The time evolution of F induced by Eq.(0) is dictated by the nonlinear 
evolution operator 
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d t r[<f>,ir;t] = -(£ d + £ q )r[0,7r;t], 



(6) 



where (fa = (0 a ,7r a )) 



cl 



and 



with 



2 2 

m 



6 



8n a (x) 



2iV 



6 (a;)0;,(x)0 a (2;) + <j> a (x)G}g(x, x) + 2<j) b (x)G^(x, x) 
— J d D xid D x 2 d D xsGat(x,xi)Gfj'(x 1 X2) 

x Gtt(x,x 3 )- 



S 3 T 



5i> i (x 1 )Silj j (x2)Sip k (x 3 ) 



S 



5ir a (x) 



8N J 1 5n b (x) 5n b (x) 5n a (x) 



G- 1 



(7) 



(8) 



(9) 



and G of (P) corresponding to G 1 ^. 

The exact flow equations for the two-point functions follow from taking 
the second derivatives of Eq.(|]) with respect to (p and n at (f> = n = 0: 



A(q) = 2uj\q)C(q) 

B(q) = -2C(q)-^j(q)B(q) (10) 
C(q) = -A(q)+^(q)B(q)-^-C(q), 



12 



where 



~ A( ^ 2) j qi q2 G(qi)G(q 2 )G(q - q 2 - Ql ) ■ 

[Au(q h -q, g 2 ) - c{q x )c{q 2 )c{q - q x - q 2 )y(q - q x - q 2 , qi, &) 

-c(gi) [2v(-g x , -g 2 , g) + v(-q x , q, -g 2 )] 

+2c(q 2 )c(-q 1 - q 2 + q)w(q - q x - q 2 , q 2 , qi)] 



= \+ ] I G{q 1 )G{q 2 )G{q-q 2 -q 1 )- 

OiV J 9i ,92 

-92) - 4c(gi)c(g 2 )c(g - gi - q 2 )z(-qx, g, -g 2 ) 
-2c(gi)w(-gi,g, -g 2 ) 

+c(g 2 )c(g -q x - g 2 ) [y(g - gi - g 2 , g 2 , -g) + 2j/(-g, q x , q-qi- ga)]] 



Similarly, the flow equation for the quartic coupling u reads []: 



^ 2 (9iM9i, 92, 93) + 4AC(gi) - 4AC(g 2 )(Si(gi + g 2 , g 3 ) 



+^ 2 (g 2 + g 3 , 9i)) - A^ 2 C(g 1 )C(g 2 )C(g 3 ) 



SYM ' 



(12) 



2 We display here only one of the six flow equations for the 4-point couplings. The 
remaining five equations can be found in Appendix A. 
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where the subscript SYM implies symmetrization with respect to the ap- 
propriate permutations of gi, q 2 , g3 and q<± = —(qi + q 2 + ^3). Here we have 
introduced the momentum integrals 

Si(?i,«0 = ^ J q G(q)G( q + qi ) ■ 

[(N + 2)u(q + q u -q, q 2 ) + 2u(q 2 , -q, q + g x ) 

-\c{q){{N + 2)v(-q, q + q u q 2 ) + 2v(-q, q 2 , q + qi)) 

+ \ c ( q + qi)c(q)((N + 2)(w(-q, q + q u q 2 ) - ^s(-q, q + q u q 2 ) 

~2 S ( q + ~ q > q2 ^ + s (~ q > q + 9l ' 92 )) 
$2(91,92) = ^ ^G(q)G(-q-q 1 )[Au(-q,q 2 ,q + q 1 ) 
+c(q)c(-q - gi)s(g + gi, -g, -gi - g 2 ) 

-c(g)(u(-g, g 2 , g + 9i) + ~9i - 92, 9 + 9i))] • (13) 

The flow equations for the quartic couplings u, v, etc. are not exact 
since we have truncated the contributions from 1PI 6-point functions. We 
furthermore have omitted the two-loop contribution to the evolution of the 
quartic couplings. Our approximation may be viewed as the second order in a 
weighted loop expansion where the evolution of every 1PI 2m-point function 
is computed in (ul + 1 — m)-loop order (i.e. two loops for the two-point 
function, one loop for the four-point function). It is easy to convince oneself 
that this expansion retains systematically all contribution in order TV 1 "™ 1 - . In 
our case it also includes (incompletely) terms of order 1/N 2 . For comparison 
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we employ a second systematic expansion, namely the 1/N expansion, where 
all terms of order 1/N 2 are omitted [] in the flow equations. 

The flow equations conserve exactly the energy density e — E/l: 



e = y^K 1 ^) + G ^)[^ + m2 + c2 ^) + ^ A / G (p)]} ( 14 ) 

N + 2 r 

^TT X / G(q 1 )G(q 2 )G(q 3 )G(q i ) [u(q u q 2 , q 3 ) ~ v(q u q 2 , gs)c(gi) 

ovv -'51,92,93 

+w(q u q 2 , q 3 )c(qi)c(q 2 ) - y(q 1} q 2 , q 3 )c(qi)c(q 2 )c(q 3 ) 
+z(qi, q 2 , q3)c(qi)c(q 2 )c(q 3 )c(q 4 )] , 

whereas the squared O(N) "angular momentum" density 

^ = N(N-1) [ G( qi )B-\ qi )- (15) 

l Jqi 

i 1 ~ In I G ^ B ~ 1 ^ q2 ^ 2w ^ q2 > ~ q2 > ~ 2s ( 9l > ~ q2 > ~ s ( ?1 > 

is conserved only up to relative corrections of order 1/N 2 . Additional inde- 
pendent conserved quantities of the form < E r (L 2 ) s > — < E > r < (L 2 ) > s 
are suppressed by inverse powers of N. They are not conserved by the trun- 
cated equations. The kinetic energy density = y J q {B~ 1 (q) + G(q)c 2 (q)} 
is, of course, not separately conserved. 

We have solved the classical flow equations (h = 0) numerically for a 
discretized system with iVj points and an ultraviolet cutoff A = a few times m, 

using a standard fourth-order Runge-Kutta algorithm which has the property 

3 With the exception of the subleading contributions to the evolution of the 4-point 
couplings that are contained in u> 2 and 7. These have to be retained in order to ensure 
exact energy conservation. 
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of being exactly time-reversible. We only consider here positive m 2 and set 
the mass scale by m = 1. For a typical cutoff A = 5 and Ni = 32, the length 
of the chain is I = N h a = vrA^/A ~ 20.1. 

3 Equilibrium properties 

As a first step in the numerical analysis, we compute the classical thermal 
equilibrium configuration (defined by the conditions: C = v = w = s = y = 
z = 0, B = (3) for different values of the parameters. In our approximation 
this corresponds to the solution of the Schwinger-Dyson equations for A and 
u that follow from the requirement d t T = 0. For X/N <C 1 it is possible to 
derive the thermal values of A and u iteratively as power series in X/N. In 
general, however, this method fails, and we found it simpler to use the flow 
equations themselves, starting from the N — > oo thermal fixed point, letting 
the system evolve for a while (At > m _1 ), taking time averages of the corre- 
lation functions, adjusting them in accordance with the thermal fixed point 
constraints, and repeating the procedure until stationary behavior with the 
desired accuracy was obtained. We were thus able to obtain configurations 
that were thermal to a very good approximation (AB/B < 0.001). In Fig. |3] 
we display the energy density and the squared angular momentum as func- 
tions of the temperature. We also show the frequency u eq (0) which is related 
in equilibrium to a (partially) renormalized temperature-dependent mass by 
<(<?) = A eq {q)T = T/G eq {q) = m\ + Z(q)q 2 . 
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1 2 3 _ 4 5 6 7 




Figure 3: Thermodynamic equilibrium properties. We show e/NT, 
L 2 /(IN(N - 1)T 2 ), and u 2 (0) as functions of T, for a system with N = 20, 
A = 2, A = 5, and I = 20.1. 
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4 Dephasing and scattering 



Next we discuss the evolution of ensembles that are initially not in thermal 
equilibrium. We opted for gaussian perturbations from the thermal state 
with initial two-point functions: 



Bo(q) 
CM 



C A A eq {q) + D, 



(9 + 9A) 



Cb\P + D b 



7? — 



(g+gBT 



D 



c 



+ e 



{q+qcr 



(16) 



The constants Da, D b , D c , qA, Qb, 1c, ^a, A b , A c are arbitrary, whereas 
Ca and Cb are tuned so that the perturbed system has the same E and L 2 as 
the unperturbed thermal equilibrium ensemble. We also use superpositions 
of gaussian perturbations with the property that the initial deviations Da, 
Db for small and large q 2 are small. 

In order to assess the importance of "scattering" for the equilibration of 
different physical quantities, we first compare the results obtained by using 
the full equations (ITUp-(O) with those obtained by keeping only the leading 
terms in 1/N (i.e. neglecting all 4-point functions). We repeat here that in 
leading order 1/N scattering is absent and only kinetic dephasing can induce 
a smoothening and averaging out of the perturbation. 

As we can see in Fig.|I|, even in the absence of interactions, energy equipar- 
tition is achieved to a very good approximation. Also u;(0) equilibrates ap- 
proximately (Fig.H). The individual correlation functions A(q) and B(q), 



however, do not equilibrate in the absence of scattering (FigsfJ, 0). They 
oscillate around constant values. When the effect of the time-dependent 
four-point functions is added the picture changes. We now see that the orig- 
inal perturbations in the correlation functions are damped and smoothened 
out by the evolution, although they do not reach exact thermal equilibrium. 
In Fig.^| we show the time evolution of A(0) and A(A/2). We remind that 
\JTA(0) should approach uj(0) in thermal equilibrium. 

5 Asymptotic behavior and thermalization 

We have seen in the previous section that next-to-leading order terms in 
the weighted loop or 1/N expansions induce an energy exchange between 
different Fourier modes, which is a prerequisite for thermalization. Also 
"particle numbers" for individual Fourier modes are no longer conserved 
separately. Because of this energy exchange, a system with nonthermal initial 
conditions is driven towards thermal equilibrium. At late times its correlation 
functions oscillate around mean values that are "more thermal" than in the 
case of mere dephasing. In Fig.[j] we show the evolution of the time-averaged 
correlation functions A(q) and B(q) in a typical case. One clearly observes 
the initial approach towards the equilibrium values. For large t, however, 
stationary values are reached which deviate from the equilibrium correlations. 
These stationary values correspond to exact fixed points of the truncated 
evolution equations. We have computed the fixed points by methods similar 
to the computation of the equilibrium. They are also displayed in Fig.[7[ 
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Figure 4: Evolution of the frequency w 2 (0), in the leading (upper) and next- 
to- leading 1/N approximations (same parameters as in Figs. 1-2). Compar- 
ison with the equilibrium value (also shown) indicates the "more thermal" 
behavior due to the inclusion of scattering, see also Fig.^j for T = 2. 
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Figure 5: A comparison between the leading-order (left) and next-to-leading- 
order (right) evolutions of the two-point correlation functions. The first row 
shows A(q)/(Yjq' B(q')/N{) — (q 2 + m 2 ), which is a measure for the deviation 
of the inverse propagator from the classical value. The second row gives 
B(q), or the deviation from the Maxwell velocity distribution B(q) = (3. The 
correlation functions are averaged over various time intervals. The initial 
values at t — are also shown. 



21 




Figure 6: Full time evolution of the two-point functions A(0) (both lead- 
ing and next-to-leading order) and A(A/2) (next-to-leading only), without 
averaging. 
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We notice that the effect of the deviation from thermal equilibrium on the 
correlation function in coordinate space is small, as can be seen from the plot 
of G(x — y). 

In Table 1 we collect the asymptotic stationary values for several differ- 
ent choices of physical parameters and temperature, and otherwise identical 
initial conditions^. They clearly differ from thermal equilibrium, typically on 
the 10% level for not too small values of X/N. As a general rule, the larger 
X/N, the faster the system approaches its asymptotic limit, and the closer 
this limit is to the thermal values. In order to save on computer time, we 
have therefore opted for rather small values of iV (1 < N < 5) and values of A 
between 1 and 60. Of course, for small iV and/or large X/N the applicability 
of a 1/N or weigthed-loop expansion is questionable. There is no satisfac- 
tory way to assess the reliability of either truncation in the strong coupling 
regime. However, we feel comforted by the fact that the 4-point functions 
never grow so large as to make the system unstable, and that the time fluc- 
tuations of L 2 (which we remind is conserved exactly by the exact evolution 
but only up to relative corrections ~ 1/N 2 by the truncated equations) are 
always quite small (at the 1% level for A = 60, N = 1). The most direct test 
seems to be a comparison between the results from the weighted-loop and 
the 1/N expansions. For N — 1, large A (Table lb), the two truncations lead 

to very different large time stationary values, and it is therefore conceivable 

This holds up to discretization corrections, since we define identical initial conditions 
by identical continuous functions of q. 
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that the asymptotic departure from thermal equilibrium is due to truncation 
errors. As N gets larger and N smaller, however, the two expansions agree 
much better (Fig||) and it is less plausible that higher order 1/N corrections 
could account for the asymptotic departures from equilibrium. We there- 
fore believe that our numerical results support the existence of nonthermal 
attractive fixed points also for the exact (i.e. non-truncated) system. The 
implications of this claim are discussed in the next section. 

As a final word of caution, we should consider the possibility that full 
thermalization does occur even for small X/N, but on longer time scales than 
those we have been able to probe. This issue can be settled only with further, 
more computer-intensive investigations. What we have seen so far does not 
support this hypothesis. Even if that turned out to be the case, our method 
would still be very useful for identifying the various time scales (dephasing - 
partial thermalization - complete thermalization) and for assessing the role 
of thermalization in practical problems, where extremely large time scales 
are not always relevant. 

We have also studied the volume dependence of the large time behavior 
for a particular choice of parameters (Table lc, Fig.^J). Apparently, the 
large volume limit still differs from thermal equilibrium. Relatively small 
volumes (I ~ 20m -1 ) seem often to be sufficient for an extrapolation to 
infinite volume. 

Finally, we have investigated a limited sample of initial conditions that 
are not symmetric under time reversal, and whose backward and forward 
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evolutions therefore differ in their microscopic details. In all cases, the large 
time asymptotic averaged values of the correlation functions come out the 
same in both time directions. 



B(A/2)-/3 A(A/2)-A e ,(A/2) 
P A eq (A/2) 



(*) 3 


1 


0.3 


20.1 


0.083±0.005 


0.08 


(*)1 


2 


0.3 


20.1 


0.158±0.005 


0.16 


(*)1 


5 


0.3 


20.1 


0.16±0.01 


0.16 


10 


1 


0.3 


20.1 


0.02±0.005 


0.03 


;*) 10 


1 


0.3 


20.1 


0.03±0.005 


0.03 


60 


3 


0.3 


20.1 


0.07±0.01 


0.07 


60 


1 


0.3 


20.1 


-0.013± 0.003 


0.013 


;*) 60 


1 


0.3 


20.1 


0.008± 0.003 


0.009 


(*)1 


3 


0.3 


5.03 


0.23±0.003 


0.23 


(*)1 


3 


0.3 


10.05 


0.20± 0.015 


0.20 


(*)1 


3 


0.3 


15.08 


0.183± 0.015 


0.24 


(*)1 


3 


0.3 


17.6 


0.18 ± 0.005 


0.18 


1 


3 


0.3 


17.6 


0.187± 0.006 


0.185 


(*)1 


3 


0.3 


20.1 


0.11± 0.03 


0.12 


1 


3 


0.3 


20.1 


0.167± 0.008 


0.17 


(*)1 


3 


0.3 


25.1 


0.167±0.017 


0.16 



Table 1: Asymptotic displacement from thermal equilibrium for different A, N, 
(3 and I. For all configurations A = 5. The initial perturbation is a superposition 
of three gaussians with D B = (3/2, q B = A/2, A B = 1.5A/16; D B = -(3 /4, 
q B = 5A/16, A B = 1.5A/16; and D B = -(3/4, q B = 11A/16, A B = 1.5A/16. 
Configurations marked by (*) are evolved according to the 1/N expansion, the 
others according to the weighted loop expansion. For the first entries (a) the 
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1/N expansion seems reliable. The second group of entries (b) concerns large 
interactions with a rapid approach to asymptotic behavior. Finally, the last entries 
(c) are used for a study of the volume dependence, see also Fig.^. 

6 Discussion 

Our study of the evolution equations, applied to various initial nonthermal 
probability distributions, clearly establishes effectively irreversible behavior. 
This is not put in by hand in the form of irreversible evolution equations. 
Our equations are manifestly invariant under time reflection. Effective irre- 
versibility is rather related to the existence of stationary solutions or fixed 
points towards which the flow is effectively attracted. It can be observed by 
evolving in time directions t — ► ±00. We see both the effects of dephasing, 
i.e. effective loss of phase information, and scattering, i.e. energy exchange 
between different momentum modes. Our investigations are carried out for 
translation invariant ensembles such that the energy exchange is not merely 
due to a classical background field evolving in time. The inclusion of scat- 
tering effects is a crucial step beyond the leading 1/N- approximation used 
in the past 0. Genuinely, the system approaches asymptotically for large t 
an oscillatory behavior of the correlation functions around a stationary so- 
lution. The time averaged values of the correlation functions are close to 
the corresponding stationary solutions. In a rough sense, the stationary so- 
lutions share many properties of thermal equilibrium. The corresponding 
fixed points are, nevertheless, not identical to the thermal fixed point. The 
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Figure 7: Time evolution of the two-point functions for a system with A = 10, 
N = 1, (3 = 0.3, with nonthermal initial conditions. The three panels show, 
from top to bottom, the relative deviations of A(q), B(q) and G(x), averaged 
over various time intervals, from their respective thermal values. 
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Figure 8: Time evolution of B(A/2) averaged over At = 30. The parameters 
are A = 1, N — 3, (5 — 0.3, I = 17.6, and the initial perturbation is described 
in the the caption of Table 1. The two curves correspond to the 1/N and 
weighted loop expansions (bold and plain, respectively). There is no sign or 
an asymptotic approach to the precise thermal value B(A/2) = 0.3. 
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Figure 9: Volume dependence. We show the asymptotic time-averaged value 
of 5 (A/2) as a function of I, for the system of Table lc. Diamonds and 
crosses correspond to the 1/N and weighted loop expansions, respectively. 
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latter turns out to be a point in a whole manifold of fixed points. For generic 
fixed points in this manifold the correlation functions differ from thermal 
equilibrium. 

Perhaps the most interesting observation concerns the difference of the 
asymptotic stationary ensembles from the thermal ensemble. The system 
retains memory of the initial conditions beyond the energy density or tem- 
perature. Isolated systems seem to differ in this respect from systems coupled 
to a heat bath. Arguments why isolated systems do not thermalize exactly 
can be given on different levels. First there are exact obstructions from con- 
served correlation functions. An example is the squared angular momentum 
density < L 2 > /I. In thermal equilibrium this quantity can be computed as 
a function of temperature, L 2 eq /l = N(N — 1)T J q G eg (q;T). Since < L 2 > is 
conserved by the exact flow equations, any initial value of < L 2 > different 
from the thermal one implies immediately that the correlation functions ap- 
pearing in Eqs.(|T^)-([T5|) cannot all take thermal values for t —>■ oo, not even 
in a time-averaged sense. We emphasize that this obstruction is based on an 
exactly conserved quantity and therefore cannot be an artifact of insufficient 
appr oximat ions . 

In principle, one could take care of the conservation of < L 2 > by an 
extension of the thermodynamic description, adding a chemical potential for 
L 2 . The problem is, however, that there exist infinitely many conserved 
combinations of correlation functions. Another prominent example is the 
global "specific heat" cy =< (E— < E >) 2 > / T 2 l which corresponds 
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again to an exactly conserved combination of correlation functions < (E— < 
E >) 2 >= 4r J q Gl n (q) + .... Also, the additional "chemical potentials" 
would multiply nonlocal expressions like E 2 = J dxdy e(x)e(y). Indeed, the 
probability distribution p ~ e -K E / l ) E j s stationary for an arbitrary function 
6(e). The Boltzmann distribution 6(e) = f3 is only a special case. If the 
function 6(e) is not constant the value of cy typically differs from the one 
in thermal equilibriumF]. These general considerations hold for an arbitrary 
number of space dimensions. They show that for isolated systems there 
cannot be a proof of strict thermalizationP] using arguments of ergodicity. 
Strict thermalization for arbitrary initial conditions is in contradiction with 
the existence of conserved combinations of correlation functions]]. 

On a second level one observes the existence of a large manifold of fixed 
points or stationary solutions for the truncated flow equations |J. In our 
truncation they are given byC = v = y = and by solving the remaining 
equations d t C = d t v = d t y = 0. The latter equations determine the station- 
ary values for A, B, u, w, s, z only incompletely. The present work clearly 
establishes numerically that these fixed points differ, in general, from the one 

corresponding to thermal equilibrium. They also prove stable with respect to 

5 The infinite volume behavior (< E 2 > — < E > 2 )/ < E 2 >~ Z -1 holds for a wide 
class of 6(e) if its deviation from a constant scales properly with I. 

6 By "strict thermalization" we mean an asymptotic approach of the probability distri- 
bution to the Boltzmann distribution, the distribution of the microcanonical ensemble. 

7 The general problem with ergodicity arguments is that only a finite neighborhood of a 
given point in phase space will be reached by an arbitrary trajectory after a finite lapse of 
time. This is not enough since even very close trajectories typically separate substantially 
at later time. 
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small fluctuations. By investigating initial conditions with the same < L 2 > 
as for thermal equilibrium we also establish explicitly that the fixed points 
are not not fully specified by < E > and < L 2 >. Numerically, we actually 
find a large manifold of different fixed points for given < E > and < L 2 >, 
as suggested by the counting of equations and variables for the general sta- 
tionary solutions. We also establish that the nonthermal fixed points play 
a role in the asymptotic dynamics. Nonthermal initial conditions typically 
result in fluctuations around nonthermal stationary solutions at late time. 

One may ask if the existence of these fixed points could not be an artifact 
of the truncation. Three arguments indicate that this is not the case. First, 
some coordinates in the fixed point manifold are related to exactly conserved 
combinations of correlation functions like < L 2 > /I. Second, the counting 
of equations and variables indicates that the dimension of the fixed point 
manifold further increases once 1PI six-point functions or higher couplings 
are included. Third, we find similar fixed points for different truncations 
in next-to- leading order in 1/N. The ones approached by a given initial 
condition are close to each other for small X/N and stay substantially away 
from the thermal fixed points. 

Our investigation of the volume dependence indicates that the fixed point 
manifold does not shrink to the thermal fixed point in the infinite volume 
limit. Furthermore, all numerical results suggest that the nontrivial fixed 
points play indeed a dynamical role. From all this a picture for the asymp- 
totic late time behavior of isolated systems emerges where some features of 
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" mesoscopic" dynamics survive even in the infinite volume limit. The initial 
information is not lost completely as for a thermalizing system. Part of the 
information survives and specifies the stationary solution around which the 
system oscillates asymptotically. 

We observe that for large enough interactions the deviations from thermal 
equilibrium are small - typically a few percent for the correlation functions 
and even less for quantities which involve momentum averages. Part of this 
can be explained by exact relations which hold for all stationary solutions. 
As an example, let us consider the condition for a static < qp > correlation 

~d t J dx < q a {x)p a (x) >=^Jdx< p a (x)p a (x) > (17) 
1 ( 

-- J dx {< diq a (x)diq a (x) > +m 2 < q a (x)q a (x) > (18) 
+ 7^ < (q a (x)q a (x)) 2 >| =0 

which relates the kinetic and potential energy (note that they are not equal 
for interacting systems): 

A c 

< E kin >=< E pot > +— J dx< (q a (x)q a (x)) 2 > . (19) 

This is, of course, the thermal relation, but it extends to all other station- 
ary solutions as well. We conjecture that large interacting systems generi- 
cally show an effective irreversible evolution towards asymptotic oscillations 



around one of the stationary solutions. Then relations of the type (|T^) hold 
asymptotically irrespectively of the initial conditions. This explains the ro- 
bustness of a large set of asymptotic time averages of correlation functions 
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- an important part of the initial information is indeed lost. Conversely, 
a judgement of precise asymptotic thermalization should not be based on 
generic relations like Eq . fll9|) , but rather on correlation functions which can 
differ for two inequivalent fixed points. 

The lack of exact thermalization of large interacting systems has conse- 
quences for "systems in a heat bath" as well. Indeed, we may consider a 
subsystem, say q a {x), p a {x) for \x\ < Iq/2 <C 1/2 and view it as evolving in 
the "heat bath" consisting of the degrees of freedom with Iq/2 < \x\ < 1/2. 
Can we expect that the subsystem effectively thermalizes even though the 
large isolated system (subsystem and bath) does not exactly thermalize? 
This question can be addressed by an investigation of correlation functions 
for the subsystem, say < q a (x)q a (y) > c with \x\, \y\ < Iq/2, or a convenient 
smoothened version (&o = tt/Iq) 



G ko (x,y) = i < q k a °(x)q k a °(y) > c = G(x - y )e-^^\ (20) 



9a W = Qa(x)e 



(21) 



For a translationally invariant ensemble one has 




(22) 




D 
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In particular, the Fourier transform of Gk (x, 0) reads 




(q-p) 2 



GkM 



G fc0 (g,g / ) = (2vr) 



G(p)e 



(23) 



Only for k$ — > does this reduce to G(q). For nonzero ko (corresponding to 
a subsystem), however, Gk {q) involves a momentum averaging with width 
given by k$. Because of dephasing, the momentum averaged two-point func- 
tion Gk {q) approaches a stationary value much more efficiently than G(q). 
For large lo, and nevertheless lo <C /, it is conceivable that Gk (q) actually 
reaches asymptotically a stationary value whereas G(q) fluctuates around a 
stationary value. Nevertheless, the time averaged values of Gk {q) and G(q) 
are still related by Eq . (|23f) . A nonthermal asymptotic behavior of the time 
averaged G(q) will manifest itself also in the asymptotic form of Gk {q) if it 
extends over a momentum range with width larger than k . Only variations 
of G(q) in small momentum ranges will be washed out. From our present in- 
vestigation we see no indication that asymptotically Gk (<?) reaches precisely 
its thermal value. We conclude that mesoscopic dynamics may also be of rel- 
evance for subsystems which are in thermal contact with a "heat bath" . The 
crucial point here is that the heat bath itself is not precisely thermalizing. 

In summary, our investigation indicates that isolated systems roughly 
thermalize for large time, while some quantitative deviations from thermal 
equilibrium remain. The "loss of memory of the initial conditions", usually 
assumed in the picture of thermalization, turns out not to be complete. This 
holds for interacting systems and in the large volume limit. Our results 
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question Boltzmann's thermalization conjecture for isolated systems. They 
suggest that even large interacting isolated systems do not thermalize in a 
strict sense. 

Our treatment is based on an exact evolution equation for the time de- 
pendence of equal time correlation functions. Nevertheless, the solution of 
these equations involves approximations in the form of a truncation of the 
time dependent effective action. Since our findings touch the basics of ther- 
modynamics, they should be questioned by an independent method. One 
possibility seems the numerical solution of the microscopic equation (|2|) for 
a large sample of different initial conditions. Taking an ensemble average 
over the initial conditions gives directly the equal time correlation functions 
which can be compared with the present work. Such a computation could 
establish definitely if the findings of this work are substantially affected by 
the truncation or not. 

Acknowledgement: We thank L. Bettencourt for many helpful discus- 
sions. 

7 Appendix A: Flow Equations 

In this appendix we present the evolution equation for the 1PI 4-point func- 
tions which are not specified in the main text. 

v(qi,q2,qs) = [2uj 2 (q 2 )(w(q u q 2 , q 3 ) ~ s(q u q 2 , q 3 ) ~ s(qi, q 2 , &)) 

+2w 2 (g 3 )s(gi, q 3 , q%) - 4u(gi, q 2 , q 3 ) + AXB(q 1 ) - -^-v(q!, ?2 
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%i,?2,ga) 



y(?i,?2,?a) 



-4A( J B(gi)(S'i(gi + g 2 , 93) + S 2 (qi + 93, 92)) + C(q 4 )S 3 (qi + q 2 , qi) 
+2C(q 2 )S 5 (q 2 + q 3 , 9l ) + C(g 4 )S 6 (gi + 93, 9i)) 
-Xh 2 B( qi )(C(q 2 )C(q 3 ) + C(g 3 )C(g 4 ) + C(g 2 )C(g 4 ))] syM 

w(qi, q 2 , q 3 ) = [u 2 (q 3 )(y(qi, 92, 9s) + y(qi, 9s, 92) + 1/(92, 9s, 9i)) - 92, 9s) 

7(91) ~l~ 7(92) 

-^(92, 94, 9i) - ^(92, 93, 9i) M9i, 92, 9s) 

-A(C(g 3 )5 , 4 (9i, 92) + SB(q 2 )S 5 (q 2 + q 3 , 
-4A( J B(gi)S , 3 (-gi - g 3 , g 2 ) + B(q 1 )S 6 (q 2 + q 3 , q 2 ) 
+C(q 3 )S 7 ( qi ,q 2 )) - 3\h 2 B( qi )B(q 2 )C(q 3 ) 



SYM 

^ 2 (q 3 )y(qi, 93, 92) - 2v(q 2 , 94, 9l) ~ 7 ^ 1 ^ 7 ^ 2 ^ g(g 1; g 2? g 3 ) 
-4A( J B( gi )5 3 (-9i - 93, 92) + B{ qi )S 6 {q 2 + q 3 , q 2 ) + C{q 3 )S 7 {q u q 2 )) 
-2Xh 2 B( qi )B(q 2 )C(q 3 )] sYM 

4uj 2 (q 4 )z(q l7 q 2 , q 3 ) - 2w(q x , q 2 , q 3 ) - s(q 2 , q 3 , g 4 ) - s(q 1 , q 3 , q 2 ) 
+s(9i, 92, 9s) + «(9i, 92, 94) - 7 ^ +1 ^ y{q u 92, 9s) 

-4A(^(g 3 ) 54(g ! ,g2) + B( qi )S 7 (q 2 , q 3 )) - Xh 2 B( qi )B(q 2 )B(q 3 ) 



%i,92,9s) = 



-y(9i, 92, 9s) - 4 7 -^-z(gi, g 2 , g 3 ) 



TV 



SYM 

(24) 



J SYM 



They involve the following momentum integrals 
1 

N + 2 



S3(9i,92) = ± f G(q)G{q - <h) 



^(92, -9, 9 - 9i) + ^(92, 9-9i, -9) 



f- 2 ^(92, 9i - 92, -9) + 2 C (^) C (^ _ Qi)(( N + 2 )y(l ~ 9i, -9, 92) 
+y(92, 9 - 9i, -9) + 1/(92, -9, 9 - 9i)) - c(q)s(-q, q 2 , q - qj 
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£4(91,92) 



£5(91,92) 



Se (91,92) 



£7(91,92) 



~c(q - qi)(Ns(q - 91, 92, -9) + 2w(g - 91, 92, -9))] 
i J q G(q)G(q- 91-92) ■ 

(iV + 2)(w(9i, 92, 9 - 9i - 92) - ^s(qx, g 2 , q - qi - q 2 ) - -s(?i, 92, -9)) 
+s(gi, 92, -9) + 2c(9 - 91 - 92)0(9) (Nz(q - 91 - 92, -9, 9i) 
+2^(92, 9 - 9i - 92, 9i) + 2z(92, 9i, 9-91-92)) 
-c(g) (JVy(gi, 92, -9) + 2y(-g, 91, 92) + 2t/(gi, 92, -9))] 

[2^(92, -9, 9 + 9i) + 20(9)0(9 + 91)2/(92, ~9, 9 + 9i) 

1 1 
-ic(q)(w(-q, 92, 9 + 9i) - 2 S(-(? ' ?2 ' 9 + 9i) - g2 ' ~ ^ 

-2c(-9 - 9i)s(g 2 , 9 + 9i, -9)] 

2^ I ft) ■ 

b(92, -9 - 9i, 9) - c(g)s(g 2 , -9, -9 - 91) - 2c(-g - qi)(w(-q - q u g 2 , 9) 
-^ s (-9 - 9i, 92, 9) - 2~ s( ~9 ~ 9i> 92, 9i - 92)) 
+c(g)c(-g - 91)2/(92, -9i - 9, 9)] 

^oMOh -„-„). 

[s(92, 9i, -9) - c(q)y(-q, 92, 9i) - c(g - 9i - 92)2/(9 - ?i - 92, 9i, 92) 
+4c(g)c(g - 91 - 92)2(91, 9 - 9i - 92, 92)] • (25) 
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